Chapter 07

1 Exercise 01

Consider the the number of pigs slaughtered in Victoria, available in the aus_livestock dataset.

1. Use the ETS() function to estimate the equivalent model for simple exponential smoothing. Find the optimal values of \(\alpha\) and \(\ell_0\), and generate forecasts for the next four months.

## Series: Count 
## Model: ETS(A,N,N) 
##   Smoothing parameters:
##     alpha = 0.3221247 
## 
##   Initial states:
##      l[0]
##  100646.6
## 
##   sigma^2:  87480760
## 
##      AIC     AICc      BIC 
## 13737.10 13737.14 13750.07

2. Compute a 95% prediction interval for the first forecast using \(\hat{y} \pm 1.96s\) where \(s\) is the standard deviation of the residuals. Compare your interval with the interval produced by R.

auto_interval <- fc %>% 
    hilo() %>% 
    unpack_hilo("95%") %>% 
    select(4:7) %>% 
    slice(1)

auto_interval

## # A tsibble: 1 x 5 [1M]
##    .mean                  `80%` `95%_lower` `95%_upper`    Month
##    <dbl>                 <hilo>       <dbl>       <dbl>    <mth>
## 1 95187. [83200.06, 107173.1]80      76855.     113518. 2019 sij

Let’s get manually those values:

sd_res <- fit %>% 
    augment() %>% 
    pull(.resid) %>% 
    sd()

auto_interval$.mean[1] + c(-1, 1) * (qnorm(0.975) * sd_res) 

## [1]  76871.35 113501.77

Well, almost close …

2 Exercise 02

Write your own function to implement simple exponential smoothing. The function should take arguments y (the time series), alpha (the smoothing parameter \(\alpha\)) and level (the initial level \(\ell_0\)). It should return the forecast of the next observation in the series. Does it give the same forecast as ETS()?

I’ve got an inspiration from this source. We are both getting the same values for this sub-exercise, but for the next one, this function yields better values.

exp_smoothing_manual <- function(y, arg_alpha, arg_ell_zero, bool_forecast_h1 = FALSE) {
    
    if (bool_forecast_h1) {
        total_len <- length(y) + 1
    } else {
        total_len <- length(y)
    }
    
    
    anti_alpha <- (1 - arg_alpha)
    
    store_predictions <- rep(NA, total_len)
    store_predictions[1] <- arg_ell_zero
    
    for (i in seq_along(store_predictions)) {
        
        if (i == 1) {
            last_y <- store_predictions[i]
            last_yhat <- store_predictions[i]
        } else {
            last_y <- y[i - 1]
            last_yhat <- store_predictions[i - 1]
        }
        
        term_01 <- arg_alpha * last_y
        
        term_02 <- anti_alpha * last_yhat
        
        yhat <- term_01 + term_02
        
        store_predictions[i] <- yhat
        
        
    }
    
    out <- yhat[length(yhat)]
    
    return(out)
    
}

fit_model_pars <- coef(fit) %>% 
    select(term, estimate)

value_manual <- exp_smoothing_manual(
    y = aus_pigs$Count,
    arg_alpha = fit_model_pars$estimate[fit_model_pars$term == "alpha"],
    arg_ell_zero = fit_model_pars$estimate[fit_model_pars$term == "l[0]"],
    TRUE
)

value_auto <- fc$.mean[1]

value_manual == value_auto

## [1] TRUE

Great!

3 Exercise 03

Modify your function from the previous exercise to return the sum of squared errors rather than the forecast of the next observation. Then use the optim() function to find the optimal values of \(\alpha\) and \(\ell_0\). Do you get the same values as the ETS() function?

Note: this source got value of 0.299 for alpha and 76379.265 for level.. But, the code below is heavily modified and doesn’t get the values of the source. Anyways, correct values are:

true_values <- coef(fit) %>% 
    select(term, estimate)

true_values

## # A tibble: 2 × 2
##   term    estimate
##   <chr>      <dbl>
## 1 alpha      0.322
## 2 l[0]  100647.

Let’s try to get optimal values:

optimize_exp_smoothing <- function(pars = c(arg_alpha, arg_ell_zero), y) {
    
    out_predicted <- rep(NA, length(y))
    
    for (i in seq_along(out_predicted)) {
        
        out_predicted[i] <- exp_smoothing_manual(
            y = y[1:i],
            arg_alpha = pars[1],
            arg_ell_zero = pars[2]
        )
        
    }
    
    squared_errors <- (out_predicted - y) ^ 2
    out <- sum(squared_errors) %>% sqrt()
    return(out)
    
}

optim_pigs <- optim(
    par = c(0.5, aus_pigs$Count[1]),
    y = aus_pigs$Count,
    fn = optimize_exp_smoothing,
    method = "Nelder-Mead"
)

true_values %>% 
    mutate(my_values = optim_pigs$par,
           pct_diff = (my_values / estimate) - 1)

## # A tibble: 2 × 4
##   term    estimate my_values   pct_diff
##   <chr>      <dbl>     <dbl>      <dbl>
## 1 alpha      0.322     0.322 -0.0000323
## 2 l[0]  100647.    99223.    -0.0141

Very small differences. For \(\alpha\), I am practically getting the correct value, while for \(\ell_{0}\) (starting value), I missed the mark by 1.41%.

4 Exercise 04

Combine your previous two functions to produce a function that both finds the optimal values of \(\alpha\) and \(\ell_{0}\), and produces a forecast of the next observation in the series.

exp_smooth_combine <- function(time_series) {
    
    optim_series <- optim(
        par = c(0.5, time_series[1]),
        y = time_series,
        fn = optimize_exp_smoothing,
        method = "Nelder-Mead"
    )
    
    out <- exp_smoothing_manual(
        y = time_series,
        arg_alpha = optim_series$par[1],
        arg_ell_zero = optim_series$par[2],
        TRUE
    )
    
    return(out)
    
}

var_forecast <-  exp_smooth_combine(aus_pigs$Count)
var_forecast

## [1] 95186.57

Is this same as forecasted value from fable?

## My Value: 95,186.57 | model value: 95,186.56

5 Exercise 05

6 Exercise 06

7 Exercise 07

8 Exercise 08

9 Exercise 09

10 Exercise 10

11 Exercise 11

12 Exercise 12

13 Exercise 13

14 Exercise 14

15 Exercise 15

16 Exercise 16

17 Exercise 17